# %% [markdown]
# # A 题医疗物资仓储基地的选址与运输

# %%
import numpy as np
import pandas as pd
import scipy as sp
from geopy import distance
import matplotlib.pyplot as plt
from scipy.stats import linregress

plt.rcParams['font.sans-serif'] = ['Simhei'] # 设置中文显示
plt.rcParams['axes.unicode_minus'] = False # 设置正常显示负号

# %% [markdown]
# ## Problem 1

# %%
locations = pd.read_excel('附件1：城市名称与所在医院位置.xls')

# %% [markdown]
# ### 线性回归

# %%
def linear_function(x):
    slope, intercept, r_value, p_value, std_err = linregress(locations.经度, locations.纬度)
    return slope*x + intercept

# %%
granularity = 100
# lats = np.linspace(locations.纬度.min(), locations.纬度.max(), granularity)
# lngs = np.linspace(locations.经度.min(), locations.经度.max(), granularity)
lngs = np.linspace(locations.经度.min(), locations.经度.max(), granularity)
lats = linear_function(lngs)
loc_warehouses = pd.DataFrame([lats, lngs]).transpose() 
loc_warehouses = loc_warehouses.rename(columns={0:'lat', 1:'lng'})

# %% [markdown]
# ### 回归直线上最优解

# %%
# 计算两个经纬度之间的距离
# loc1：一个元组表示经纬度地址（纬度，经度）
# loc2: 城市信息
def dis_lng_lat(loc1, loc2):
    # print(loc1, loc2)
    return distance.geodesic(loc1, (loc2.loc['纬度'], loc2.loc['经度'])).km

# %%
# 计算所有城市到仓库之间的距离
# loc_warehouse 一个元组表示经纬度地址
def get_dis(loc_warehouse, locations=locations):
    return locations.apply(lambda x: dis_lng_lat(
        (loc_warehouse.loc['lat'], loc_warehouse.loc['lng']), x), axis=1).sum()

# %%
loc_warehouses['sum_dis'] = loc_warehouses.apply(lambda x: get_dis(x), axis=1)
center_circle = loc_warehouses[loc_warehouses['sum_dis'] == loc_warehouses['sum_dis'].min()]
center_circle

# %% [markdown]
# ### 基于回归最优解为圆心扩散寻找最优解

# %%
granularity = 100
circle_lat = np.linspace(center_circle.lat.values[0]-1, center_circle.lat.values[0]+1, granularity)
circle_lng = np.linspace(center_circle.lng.values[0]-1, center_circle.lng.values[0]+1, granularity)
circle_lat = pd.DataFrame(circle_lat)
circle_lat['key'] = 1
circle_lng = pd.DataFrame(circle_lng)
circle_lng['key'] = 1
circle = pd.merge(circle_lat, circle_lng, on='key')
circle = circle.rename(columns={'0_x':'lat', '0_y':'lng'})
circle = circle.drop(['key'],axis=1)
dis = circle.apply(lambda x: get_dis(x), axis=1)

# %%
res = circle
res['distance'] = dis
best_location = res[res['distance'] == res['distance'].min()]
best_location

# %% [markdown]
# ### End & 可视化

# %%
# 计算某个点到所有城市的距离之和
# pd.DataFrame({'lng':[104.114], 'lat':[30.695]}).apply(lambda x: 
#     get_dis(x), axis=1).values[0]

# %%
regression_lngs = np.linspace(locations.经度.min(), locations.经度.max(), granularity)
regression_lats = linear_function(regression_lngs)


# %%
ax = plt.subplot()
ax.scatter(locations.经度, locations.纬度)
ax.plot(regression_lngs, regression_lats, color='y')
# ax.scatter(104.114, 30.695, c='r')
ax.scatter(best_location.lng, best_location.lat, c='r', marker='*', s=100)

ax.spines['right'].set_color('none') 
ax.spines['top'].set_color('none')
ax.spines['left'].set_color('none') 
ax.spines['bottom'].set_color('none')
ax.yaxis.set_tick_params(which='both', length=0)
ax.xaxis.set_tick_params(which='both', length=0)
ax.grid(axis='y', alpha=0.4) 
ax.grid(axis='x', alpha=0.4) 
ax.set_title('City location & Warehouses location')
ax.set_ylabel("Latitude")
ax.set_xlabel("Longitude")
ax.legend(['线性拟合', '城市', '选址'])

for city in range(locations.shape[0]-1):
    loc_city = locations.loc[city].values.tolist()
    ax.text(loc_city[1], loc_city[2], loc_city[0], ha='center', va='bottom')

plt.savefig('City location & warehouses location & line.jpg', dpi=200, bbox_inches='tight', pad_inches=0.1)
plt.show()

# %% [markdown]
# ## Problem 2

# %% [markdown]
# ### 计算城市之间的距离

# %%
locations = pd.read_excel('附件1：城市名称与所在医院位置.xls')

# %%
site_selection = pd.DataFrame(best_location).drop('distance', axis=1).rename(columns={'lat':'纬度', 'lng':'经度'})
site_selection['城市名称'] = ['Start']
locations = locations.append(site_selection)

# %%
# 计算指定城市到其它所有城市的距离
def all_dis(departure, locations=locations):
    return locations.apply(lambda x: dis_lng_lat(
        departure, x), axis=1)
    

# %%
cities = locations.set_index('城市名称')
city_dis = cities.apply(lambda x: all_dis((x[['纬度']].loc['纬度'],
    x[['经度']].loc['经度'])), axis=1)
city_name = cities.index.values.tolist()
cnt = 0
for i in range(len(city_name)):
    city_dis = city_dis.rename(columns={city_dis.columns[i]:city_name[i]})
city_dis.to_csv('城市之间的距离.csv', encoding='utf_8_sig')
city_dis.shape

# %% [markdown]
# ### 贪心求解最优路线

# %%
# Mi-26 运输直升机最大载重
max_weight = 12000
# Mi-26 运输直升机最大航程
max_distance = 2000
# Mi-26 运输直升机飞行速度
velocity = 255
# Mi-26 运输直升机数量
mi26_num = 10

# 标记该城市是否已经被送物资
vis_goods = pd.DataFrame(locations['城市名称'])
vis_goods = vis_goods.set_index('城市名称')
vis_goods['vis'] = False
vis_goods.at['Start', 'vis'] = True # 起点不用送
# 初始化距离矩阵，将0填充为inf
city_dis2 = city_dis.replace(0, np.inf)

# 模拟飞机飞行
start_loc, flag =  'Start', False
total_mileage = 0
routes = []
for mi26 in range(1, mi26_num+1):    
    now_loc = start_loc
    now_mileage = 0
    route = []
    cnt = 0
    print('---------'*8)
    while True:
        print('城市', cnt+1, '当前：', now_loc, '返航里程', city_dis[now_loc].loc['Start'], 
                '剩余里程', max_distance - now_mileage)
        # 判断飞机里程是否上限，当(返航里程+50)小于(剩余里程)的时候返航
        if city_dis[now_loc].loc['Start'] + 50 > max_distance - now_mileage or len(vis_goods['vis'].drop_duplicates()) == 1:
            print('- 当前城市', now_loc, '返航距离', city_dis[now_loc].loc['Start'])
            print('- 当前飞机里程不足，准备返航')
            print('- 本次飞行路线', route)
            print('- 此路线总里程  ', round(now_mileage, 2) ,'+', 
                round(city_dis[now_loc].loc['Start']), '=', 
                round(now_mileage + city_dis[now_loc].loc['Start'], 2), 'km')
            total_mileage += now_mileage
            routes.append(['Start'] + route)
            if len(vis_goods['vis'].drop_duplicates()) == 1:
                flag = True
                print('---------'*8)
                print("All goods have arrived")
                print('Total mileage:', round(total_mileage, 2), 'km')
            break

        now_city_name = now_loc
        # 筛选出当前城市与其它城市之间的距离
        dis = pd.DataFrame(city_dis2[now_city_name])
        # 获取离当前城市最小距离的城市
        while True:
            next_city = dis[dis[now_city_name] == dis[now_city_name].min()]
            # 判断城市是否被访问过
            if not vis_goods.loc[next_city.index[0]][0]:
                now_mileage += dis[now_city_name].min()
                route.append(next_city.index[0])
                # print('已找到目标城市  ', next_city.index[0])
                # print('    距离当前城市', dis[now_city_name].min())
                # print('    飞机已行驶  ', now_mileage)
                break
            else:
                dis = dis.drop(next_city.index[0])
        # 将这个城市标记为已经访问
        vis_goods.at[next_city.index[0], 'vis'] = True
        # mi26到下一个城市
        now_loc = next_city.index[0]
        cnt += 1
    if flag:
        break


# %% [markdown]
# ### 可视化

# %%
ax = plt.subplot()
ax.scatter(locations.经度, locations.纬度) # 绘制城市
# ax.plot(regression_lngs, regression_lats, color='y') # 绘制回归曲线
# ax.scatter(104.114, 30.695, c='r') 
ax.scatter(best_location.lng, best_location.lat, c='r', marker='*', s=100) # 绘制医疗药品仓库选址

# 美化表格
ax.spines['right'].set_color('none') 
ax.spines['top'].set_color('none')
ax.spines['left'].set_color('none') 
ax.spines['bottom'].set_color('none')
ax.yaxis.set_tick_params(which='both', length=0)
ax.xaxis.set_tick_params(which='both', length=0)
ax.grid(axis='y', alpha=0.4) 
ax.grid(axis='x', alpha=0.4) 
ax.set_title('City location & Warehouses location')
ax.set_ylabel("Latitude")
ax.set_xlabel("Longitude")

# 添加城市名称
for city in range(locations.shape[0]-1):
    loc_city = locations.loc[city].values.tolist()
    ax.text(loc_city[1], loc_city[2], loc_city[0], ha='center', va='bottom')

# 运输机飞行路线图
for route in routes:
    draw_lng, draw_lat = [], []
    for r in route:
        draw_lat.append(locations[locations['城市名称'] == r].经度.values[0])
        draw_lng.append(locations[locations['城市名称'] == r].纬度.values[0])
    ax.plot(draw_lat, draw_lng)

ax.legend(['路线一', '路线二', '城市', '选址'])
plt.savefig('City location and warehouses location.jpg', dpi=200, bbox_inches='tight', pad_inches=0.1)
plt.show()

# %% [markdown]
# ## Problem 3

# %%
loss_rate = pd.read_excel('附件3：医疗物资的损耗率.xls')
loss_rate = loss_rate.rename(columns={'医疗物资': 'time_lossRate'})
loss_rate = loss_rate.set_index('时间/分钟').to_dict()
loss_rate = loss_rate['time_lossRate']

# %% [markdown]
# ### 综合考虑损耗率和距离的计算

# %%
# Mi-26 运输直升机最大载重
max_weight = 12000
# Mi-26 运输直升机最大航程
max_distance = 2000
# Mi-26 运输直升机飞行速度
velocity = 255
# Mi-26 运输直升机数量
mi26_num = 10
# 损耗率
low_loss_rate = 0.3
# 停靠时间
stop_time = 5


# %%
# 标记该城市是否已经被送物资
vis_goods = pd.DataFrame(locations['城市名称'])
vis_goods = vis_goods.set_index('城市名称')
vis_goods['vis'] = False
vis_goods.at['Start', 'vis'] = True # 起点不用送
# 初始化距离矩阵，将0填充为inf
city_dis2 = city_dis.replace(0, np.inf)

# 模拟飞机飞行
start_loc, flag =  'Start', False
total_mileage = 0
routes = []
for mi26 in range(1, mi26_num+1):    
    now_loc = start_loc
    now_mileage = 0
    route = []
    cnt = 0
    print('---------'*8)
    while True:
        print('城市', cnt+1, '当前：', now_loc, end=' ')
        route_time = int(60*now_mileage/velocity/10 + cnt/10*stop_time)*10
        if route_time < 10:
            route_time = 10
        print('当前耗时 ', route_time, 'm', '损耗率', loss_rate[route_time]*100, '%')
        # 判断飞机里程是否上限，当(返航里程+50)小于(剩余里程)的时候返航
        if loss_rate[route_time]*100>=low_loss_rate or city_dis[now_loc].loc['Start'] + 50 > max_distance - now_mileage or len(vis_goods['vis'].drop_duplicates()) == 1:
            print('- 当前城市', now_loc, '返航距离', city_dis[now_loc].loc['Start'])
            print('- 当前药品损耗率临近最高损耗率，准备返航')
            print('- 本次飞行路线', route)
            print('- 此路线总里程  ', round(now_mileage, 2) ,'+', 
                round(city_dis[now_loc].loc['Start']), '=', 
                round(now_mileage + city_dis[now_loc].loc['Start'], 2), 'km')
            total_mileage += now_mileage
            routes.append(['Start'] + route)
            if len(vis_goods['vis'].drop_duplicates()) == 1:
                flag = True
                print('---------'*8)
                print("All goods have arrived")
                print('Total mileage:', round(total_mileage, 2), 'km')
            break

        now_city_name = now_loc
        # 筛选出当前城市与其它城市之间的距离
        dis = pd.DataFrame(city_dis2[now_city_name])
        # 获取离当前城市最小距离的城市
        while True:
            next_city = dis[dis[now_city_name] == dis[now_city_name].min()]
            # 判断城市是否被访问过
            if not vis_goods.loc[next_city.index[0]][0]:
                now_mileage += dis[now_city_name].min()
                route.append(next_city.index[0])
                break
            else:
                dis = dis.drop(next_city.index[0])
        # 将这个城市标记为已经访问
        vis_goods.at[next_city.index[0], 'vis'] = True
        # mi26到下一个城市
        now_loc = next_city.index[0]
        cnt += 1
    if flag:
        break


# %% [markdown]
# ### 可视化

# %%
ax = plt.subplot()
ax.scatter(locations.经度, locations.纬度) # 绘制城市

# 美化表格
ax.spines['right'].set_color('none') 
ax.spines['top'].set_color('none')
ax.spines['left'].set_color('none') 
ax.spines['bottom'].set_color('none')
ax.yaxis.set_tick_params(which='both', length=0)
ax.xaxis.set_tick_params(which='both', length=0)
ax.grid(axis='y', alpha=0.4) 
ax.grid(axis='x', alpha=0.4) 
ax.set_title('综合考虑损耗率和距离的路线方案')
ax.set_ylabel("Latitude")
ax.set_xlabel("Longitude")

# 添加城市名称
for city in range(locations.shape[0]-1):
    loc_city = locations.loc[city].values.tolist()
    ax.text(loc_city[1], loc_city[2], loc_city[0], ha='center', va='bottom')

ax.text(104.5, 27, '最高损耗率:'+str(low_loss_rate)+'%', ha='center', va='bottom')

# 运输机飞行路线图
cnt = 0
for route in routes:
    cnt += 1
    draw_lng, draw_lat = [], []
    for r in route:
        draw_lat.append(locations[locations['城市名称'] == r].经度.values[0])
        draw_lng.append(locations[locations['城市名称'] == r].纬度.values[0])
    ax.plot(draw_lat, draw_lng, label='Route '+str(cnt))

ax.scatter(best_location.lng, best_location.lat, c='r', marker='*', s=100) # 绘制医疗药品仓库选址

ax.legend()
plt.savefig('Distance and loss rate.jpg', dpi=200, bbox_inches='tight', pad_inches=0.1)
plt.show()


